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® Check for updates 


The rate of global-mean sea-level rise since 1900 has varied over time, but the 
contributing factors are still poorly understood’. Previous assessments found that 


the summed contributions of ice-mass loss, terrestrial water storage and thermal 
expansion of the ocean could not be reconciled with observed changes in 
global-mean sea level, implying that changes in sea level or some contributions to 
those changes were poorly constrained’, Recent improvements to observational 
data, our understanding of the main contributing processes to sea-level change and 
methods for estimating the individual contributions, mean another attempt 

at reconciliation is warranted. Here we present a probabilistic framework to 
reconstruct sea level since 1900 using independent observations and their inherent 
uncertainties. The sum of the contributions to sea-level change from thermal 
expansion of the ocean, ice-mass loss and changes in terrestrial water storage is 
consistent with the trends and multidecadal variability in observed sea level on both 
global and basin scales, which we reconstruct from tide-gauge records. Ice-mass 
loss—predominantly from glaciers—has caused twice as much sea-level rise since 
1900 as has thermal expansion. Mass loss from glaciers and the Greenland Ice 

Sheet explains the high rates of global sea-level rise during the 1940s, while a sharp 
increase in water impoundment by artificial reservoirs is the main cause of the 
lower-than-average rates during the 1970s. The acceleration in sea-level rise since the 
1970s is caused by the combination of thermal expansion of the ocean and increased 
ice-mass loss from Greenland. Our results reconcile the magnitude of observed 
global-mean sea-level rise since 1900 with estimates based on the underlying 
processes, implying that no additional processes are required to explain the 
observed changes in sea level since 1900. 


Global-mean sea level (GMSL) has increased by approximately 
1.5mm yr” (refs. ++”) over the twentieth century, modulated by large 
multidecadal fluctuations®. Changes in GMSL are the net result of many 
individual geophysical and climatological processes, with some of the 
largest contributions coming from ice-mass loss and thermal expansion 
of the ocean. The level of agreement between the sum of these individ- 
ual contributions and the observed changes in GMSL—often described 
as the ‘sea-level budget’—is a key indicator of our understanding of the 
drivers of sea-level rise’. Multiple studies show closure of the sea-level 
budget within their stated uncertainties since the 1960s and over the era 
of satellite altimetry since 19935", However, rates of GMSL change and 
their contributions to the budget over the entire twentieth century, and 
especially the first half of the twentieth century, have not yet been fully 
explained or attributed. Previous observation-based studies concluded 
that the GMSL budget for the whole twentieth century could not be 
closed within the estimated uncertainties’. Various explanations for 
this non-closure have been proposed, including an overestimation of 


the tide-gauge-derived rates of GMSL change” and underestimation 
of the ice-sheet contribution”, but there is no agreement yet on the 
cause of this discrepancy”. 

Over the past few years, revised estimates of the main known driving 
processes of global sea-level rise that cover the entire twentieth century 
have become available“, the spread among different estimates of 
twentieth-century glacier mass loss has been reduced", and improved 
mapping methods and correction of instrumental bias have resulted 
in higher estimates of the contribution from thermal expansion since 
the 1960s”. In parallel, estimates of twentieth-century GMSL change 
have converged to lower rates than previously estimated, as a result of 
improved reconstruction approaches, spatial-bias correction schemes, 
and the inclusion of estimates of local vertical land motion (VLM) at 
tide-gauge locations*®”°. As a result of these developments, the GMSL 
budget needs to be re-estimated, to determine whether the observed 
sea-level rise since 1900 can be reconciled with the estimated sum of 
contributing processes. 
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Fig. 1| Observed GMSL and contributing processes. a, Observed GMSL, and 
the estimated barystatic and thermosteric contributions and their sum. b, The 
barystatic contribution and its individual components. The TWS term is the 
sum of groundwater depletion, water impoundmentin artificial reservoirs and 
the natural TWS term. c, 30-year-average rates of observed GMSL change and of 


Estimating the sea-level budget 


To obtain estimates of changes in global ocean mass (barystatic 
changes), we combine estimates of mass change for glaciers!” ice 
sheets” and terrestrial water storage (TWS). For the TWS estimate, 
we consider the effects of natural TWS variability”, water impoundment 
in artificial reservoirs and groundwater depletion”’**. For 2003-2018, 
we use observations from the Gravity Recovery and Climate Experiment 
(GRACE)*’ to quantify the barystatic changes. We estimate changes in 
sea level due to global thermal expansion (thermosteric changes) from 
in situ subsurface observations” ” over the period 1957-2018, and com- 
bine these estimates with an existing thermosteric reconstruction”. To 
obtain an estimate of GMSL changes and their accompanying uncertain- 
ties, we combine tide-gauge observations with estimates of local VLM 
from permanent Global Navigation Satellites System (GNSS) stations 
and with the difference between tide-gauge and satellite-altimetry 
observations. 

Each tide-gauge and VLM record is affected by glacial isostatic 
adjustment (GIA) and by the effects of gravity, rotation and deforma- 
tion (GRD) from contemporary surface-mass redistribution due to 
changes inice mass and TWS. Owing tothe irregular spatial distribution 
of tide-gauge sites, these effects could bias reconstructed global-mean 
and basin-mean sea-level changes”. To avoid this bias, we remove the 
local sea-level and VLM imprints from GIA and contemporary GRD 
effects from each tide-gauge and VLM record before computing 
basin-mean and global-mean sea-level changes from the tide gauges’. 

We propagate the uncertainties and associated covariances in the 
sea-level observations, inthe contributing processes, and inthe GIA and 
contemporary GRD effects into the final estimates of sea-level changes 
and the contributing processes. To this end, we generate an ensemble 
of 5,000 realizations of global-mean and basin-mean sea-level changes 
and all of the contributing processes. For processes for which multiple 
estimates are available, such as GIA, we randomly select one of these 
estimates when computing each individual ensemble member. For 
processes for which an estimate of the uncertainty is available, such 
as GNSS observations, we sample the estimate assuming a Gaussian 
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GMSL change as a result of the different contributing processes. d, 30-year- 
average rates of GMSL change due to the barystatic contribution 

and its individual components. The shaded regions denote 90% confidence 
intervals. The values inaand bare relative to the 2002-2018 mean. 


distribution of the stated uncertainty about the corresponding mean. 
Then, we compute global-mean and basin-mean sea-level changes and 
the contributing processes for each ensemble member. We use the 
ensemble mean and spread to estimate all basin-mean and global-mean 
sea-level contributions and the associated confidence intervals. See 
Extended Data Fig. 1 and Methods for a detailed description of our 
approach. 


Global-mean sea level 


Our GMSL estimate (Fig. 1a) shows a trend of 1.56 + 0.33 mm yr” (90% 
confidence interval) over 1900-2018. It is also characterized by sub- 
stantial multidecadal variability, with higher rates of sea-level rise 
during the 1940s and since the 1990s, and lower rates around 1920 
and 1970. The higher rates at the turn of the millennium are in good 
agreement with independent satellite-altimetry observations™. The 
observed trend over 1900-2018 is consistent with the sum of the esti- 
mated thermal expansion and changes in ocean mass, which sum to 
1.52 + 0.33 mm yr” (90% confidence interval). This consistency holds 
not only for the trends over the full study period, but also over the 
past 50 years (Table 1), and for the pattern of multidecadal variability 
(Fig. 1c), except for the low rates of sea-level change around the 1920s 
and early 1930s. 

Thermosteric and barystatic sea-level changes show similar multidec- 
adal variability patterns to the GMSL changes, although the amplitude 
of barystatic variability is larger than that of thermosteric variabil- 
ity, and barystatic variability is the main cause of multidecadal GMSL 
variability (Fig. 1c). The barystatic variability is not dominated by a 
single process (Fig. 1d). The above-average rate of GMSL rise in the 
1940s is largely attributable to above-average contributions from 
glaciers and the Greenland Ice Sheet, whereas the high rate of barys- 
tatic sea-level rise since 2000 is attributable to both the Greenland 
and Antarctic ice sheets and to TWS. The low rates around 1970 are 
dominated by the TWS term (Fig. 1d). This negative contribution is 
caused predominantly by reservoir impoundment. Between 1900 and 
2003, 9,400 + 3,100 km? (90% confidence interval) of water has been 


Table 1| Linear trends in observed GMSL and in individual 
contributions to GMSL 


1900-2018 1957-2018 1993-2018 

(mm yr") (mm yr") (mm yr") 
Glaciers 0.70 [0.52, 0.89] 0.52 [0.36, 0.73] 0.67 [0.53 0.84] 
Greenland Ice Sheet 0.44 [0.35, 0.53] 0.30 [0.21, 0.38] 0.65 [0.57 0.74] 
Antarctic Ice Sheet 0.08 [0.00, 0.17] 0.13 [0.04, 0.22] 0.32 [0.21 0.44] 
TWS -0.21 [-0.34, -0.08] -0.14 [-0.31, 0.02] 0.31 [0.14 0.50] 
Barystatic 1.00 [0.71, 1.31] 0.80 [0.49, 113] 1.97 [1.63 2.33] 
Thermosteric 0.52 [0.34, 0.69] 0.71 [0.54, 0.88] 1.19 [0.99 1.44] 
Summed 1.52 [1.20, 1.85] 1.51 [1.18, 1.84] 3.16 [2.78 3.57] 
contributions 
Observed GMSL 1.56 [1.24, 1.89] 1.78 [1.48, 2.07] 3.35 [2.91 3.82] 
Observed GMSL 0.04 [-0.31, 0.41] 0.26 [-0.07, 0.59] 0.19 [-0.32 0.70] 


minus summed 
contributions 


Satellite altimetry 3.32 [2.87 3.79] 


The numbers in brackets indicate the 90% confidence interval. 


impounded, leading to a sea-level drop of 26 + 9 mm (90% confidence 
interval), witha peak in dam construction around the 1970s”°. The rate 
of global thermosteric sea-level rise since 2000 is significantly greater 
than at any momentin the twentieth century. However, the barystatic 
rate since 2000 is not significantly greater than the rate in the 1930s. 
The only major feature in observed GMSL that is not replicated by the 
sum of the processes is the low rate in observed sea-level change during 
the 1920s, although this low rate is found in most ocean basins and is 
also visible in other reconstructions (Extended Data Fig. 2). A possible 
explanation for this mismatch could be the low number of available 
tide-gauge records over the first few decades of data, which results in 
aless robust reconstruction (Extended Data Fig. 3) and in increasing 
unquantified uncertainties in individual budget components. 

The relative contributions of the barystatic and thermosteric com- 
ponents to GMSL vary over time. Figure 2a shows that the barystatic 
component dominates over the first half of the twentieth century, 
explaining more than 80% of total GMSL rise. The barystatic contribu- 
tionis larger than the thermosteric contribution over most of the sec- 
ond half of the century too, except during the peak of dam construction 
in the 1970s. Glaciers are the largest contributor to sea-level rise over 
most of the twentieth century, overtaken by the thermosteric contri- 
bution only after 1970. In Fig. 2b, we omit the TWS term to remove the 
direct anthropogenic contributions due to reservoir impoundment 
and groundwater depletion. Without the TWS term, the relative con- 
tribution from glaciers and ice sheets gradually decreases during the 
end of the twentieth century; however, their combined contribution 
increases again from the start of the twenty-first century. This increase 


is consistent with recent assessments of the sea-level budget over the 
satellite era”. 


Basin-mean sea level 


The global changes can be broken down into basin-mean changes 
(Fig. 3, Extended Data Table 1), each with different trends and vari- 
ability. Although salinity-induced (halosteric) changes in sea level cause 
negligible changes in GMSL*, they can be important contributors at 
the ocean-basin level. Thus, basin-mean changes in sea level due to 
changes in water density (steric changes) cannot be approximated by 
thermosteric changes alone*’. Because in situ salinity estimates before 
the 1950s are too sparse to extract basin-scale salinity changes, we can 
assess the basin-mean sea-level budget only since the 1950s. 

Over 1957-2018 and 1993-2018, the sea-level budget in each basin 
is closed within the 90% confidence intervals. The uncertainties of 
regional sea-level reconstructions vary considerably among basins. 
This is not only because of differences in tide-gauge coverage (Extended 
Data Fig. 3), but also to a large extent because of uncertainties in the 
GIA correction. In some basins, most tide-gauges are located in areas 
with large GIA uncertainties, such as the northwest Atlantic and the 
northeast Pacific coasts. On the other hand, the large uncertainties 
in the South Atlantic can be linked to the low number of tide-gauge 
records, with only a few records available before the 1960s. 

In contrast to the global-mean variability, which is dominated by 
barystatic variability, basin-mean multidecadal sea-level variability 
is dominated by steric changes. The steric trends vary considerably 
between basins: for example, since 1957, the subtropical North Atlantic 
has experienced asteric trend 2.7 + 0.4 times higher than the east Pacific. 
Ocean-mass trends in each basin are more homogeneous, except for 
the low trend inthe subpolar North Atlantic. This low trend is due tothe 
proximity to the Greenland Ice Sheet and regions of substantial glacier 
mass loss. Owing to GRD effects, oceans near areas of land-mass loss 
see below-average ocean-mass increases (Extended Data Fig. 4). This 
below-average increase is partially offset by GIA, which causes an upward 
trend inthis basin. As a result, despite the fact that the observed sea-level 
changes in the subpolar North Atlantic can be attributed to a different 
mix of processes, the resulting trend since 1900 is of similar magnitude 
to the global-mean. GIA also results in above-average sea-level trends 
in the subtropical North Atlantic; for other basins, its contribution is 
negligible compared to ocean-mass and steric contributions. In each 
basin, thetrend since 2000 is larger than the trend over the entire period. 
These high rates of GMSL change since 2000 are seen globally and are 
not driven by processes limited to a subset of ocean basins. 


Conclusions 


We reconstructed the GMSL since 1900 and compared it tothe sum of 
the contributing processes. We found that these processes explain the 
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Fig. 2 | Fraction of the 40-year-average summed rate explained by each contributor. a, Fraction with all components included. b, Fraction after omitting the 


TWS component. The shaded regions denote 90% confidence intervals. 
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Fig. 3 | Observed basin-mean sea level and contributing processes. 
a-f, Observed basin-mean sea level, and the estimated contributions and their 
sum, for the different basins (as indicated on the map). Contrary to the global 


observed twentieth-century GMSL trend and match the multidecadal 
variability pattern, except for the low rates in observed sea-level rise 
during the 1920s. Barystatic changes are the primary contributor to 
sea-level rise, with glacier mass loss being the largest component. Reser- 
voir impoundment caused a substantial, albeit temporary, slowdown of 
GMSL rise during the 1970s. The relative contributions of thermosteric 
and barystatic changes to GMSL vary with time. On basin scales, trends 
and multidecadal variability deviate from the global mean, mostly as 
aresult of variability in the steric component. 

In the subpolar North Atlantic, along which almost half of all tide 
gauges used in this study are located, including many of the longest 
available records, the ocean-mass contribution over the twentieth 
century is negligible, whereas GIA causes relative sea level to rise in 
this basin. This combination results in sea-level trends that are compa- 
rable to global-mean trends, but caused by a different combination of 
processes. Although many of the world’s longest tide-gauge records, 
including the 225-year record from Amsterdam and the 220-year record 
from Brest, are located along the coast of the subpolar North Atlantic, 
long-term changes derived from these records are not representative 
of global-mean changes. 

Closure of the twentieth-century sea-level budget, as demonstrated 
here, implies that no additional unknown processes, suchas large-scale 
deep-ocean thermal expansion or additional mass loss from the Ant- 
arctic Ice Sheet, are required to explain the observed changes in global 
sea level. Such additional processes had been speculated to explain the 
non-closure found in previous studies of global sea-level budget”*”. Our 
demonstration of closure of the global-mean and basin-mean sea-level 
budget forms a consistent baseline against which process-based and 
semi-empirical sea-level projections can be benchmarked, without 
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case, GIA causes basin-mean changes in sea level, and so is included inthe sum 
of contributors. The shaded regions denote the 90% confidence interval. The 
values are relative to the 2002-2018 mean. 


the need to compare against either the sum of processes or observed 
sea level”. The downward revision of the estimated sea-level rise 
and updated estimates of the driving processes, particularly the 
increased estimated glacier mass loss, result in a consistent picture 
of twentieth-century GMSL rise and its underlying causes. 
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Methods 


The global-mean and basin-mean sea-level changes that we report are 
relative sea-level (RSL) changes®, corresponding to the total change 
in sea-water volume. RSL changes are changes relative to the underly- 
ing seafloor. They differ from geocentric sea-level changes observed 
by satellite altimetry, owing to seafloor deformation. We divide the 
global ocean into six basins’. These basins (Extended Data Fig. 3) are 
defined using a clustering approach that merges locations that share 
a common interannual sea-level variability signal, as observed by 
satellite altimetry. We define the global ocean as the sum of all basins. 
Our basins do not cover the highest latitudes of polar oceans, as 
satellites cannot sufficiently provide data for these regions. Sea-level 
changes in these regions, which cover 7% of the total ocean area, are not 
included. Because the omitted area is small, only alarge local anomaly 
in sea-level rise would have to potential to affect GMSL substantially. 
A recent sea-level reconstruction’ estimates a rate of sea-level rise of 
1.0 + 0.8 mm yr” in the Arctic ocean and a rate of 1.6 + 0.6 mm yrtin 
the Southern Ocean over 1900-2015. Using these rates to extend our 
reconstruction has a negligible (less than 0.1mm yr”) effect on the 
global-mean sea-level trend. Therefore, omitting these oceans when 
reconstructing global-mean sea-level changes is unlikely to cause 
substantial GMSL changes. 


The ensemble approach 

Assessing closure of the global-mean and basin-mean sea-level budget 
requires an estimate of the mean and associated uncertainties of the 
observed sea-level changes, as well as those of the major contributing 
processes. Some processes, especially GIA, affect both the sea-level 
observations and estimates of the contributing processes, and the 
reconstructed sea-level changes and the sum of processes are not fully 
independent. Therefore, we use a Monte Carlo approach to obtain a 
consistent set of observed sea level, its contributing processes and 
associated uncertainties. We generate 5,000 realizations of observed 
sea level and the contributing processes. For each process, we use one 
of the two following approaches. If a large number of estimates is avail- 
able, we randomly select one estimate (for example, GIA). If only a 
single or limited number of independent estimates are available (for 
example, glacier mass loss), we generate ensemble members by ran- 
domly selecting and perturbing one of these estimates. We perturb the 
estimate by drawing random numbers from a Gaussian distribution 
using the a priori uncertainty of that estimate as the standard devia- 
tion and adding these random numbers to the estimate. We compute 
basin-mean and global-mean sea-level changes and the contributing 
processes for each ensemble member. This procedure provides 5,000 
realizations of global-mean and basin-mean sea level, allcomponents, 
and the difference between sea level and the sum of the components, 
in which all known sources of uncertainty and the spread among differ- 
ent estimates have been propagated. We compute all the time series, 
moving trends and linear trends for each ensemble member and subse- 
quently derive the mean and confidence intervals from the ensemble. 
This procedure ensures that the underlying co-variances between the 
sea-level observations and contributing processes are propagated into 
the final estimates. Extended Data Fig. 1 shows the procedure that is 
followed for each individual ensemble member. In the sections below, 
we describe the data and estimates used for reconstructing sea level 
and each process. 


GIA 

While not changing contemporary GMSL, GIA causes changes in the 
Earth’s gravity field and the shape of the solid Earth, and changes local 
relative sea level. These changes affect observations from tide gauges, 
altimetry, GNSS stations and our estimates of the contributors to bar- 
ystatic sea-level changes”. Estimates of GIA-induced changes in sea 
level, gravity and the solid Earth all come witha substantial uncertainty. 


Because GIA input parameters simultaneously affect several compo- 
nents of the sea-level budget, these components and their uncertainties 
are not fully independent of each other**°. To estimate the GIA effects 
and to propagate the mutually dependent uncertainties in the GIA 
predictions into all affected observations, we use an ensemble of 
GIA estimates“. This study provides a128,000-member ensemble of 
GIA predictions, computed by varying solid-Earth parameters (litho- 
sphere thickness and mantle viscosities) and amplitudes of global 
deglaciation histories over the past 20,000 years. Each GIA ensem- 
ble member provides a consistent set of changes in relative sea level, 
solid-Earth deformation and changes in equivalent water height, 
used to correct GRACE observations, and comes with a likelihood 
that reflects how good the fit is to a dataset of vertical GNSS veloci- 
ties and palaeo sea-level records. Therefore, this model allows for a 
robust quantification of the uncertainties associated with GIA. The 
spread between the ensemble members depicts the uncertainty inthe 
GIA predictions due to uncertainty in the solid-Earth parameters and 
the deglaciation history. Large uncertainties can therefore be found 
around the edges of formerly glaciated regions, suchas the coastlines 
of Alaska and Fennoscandia, and the forebulge collapse regions along 
the North American coastlines. The ensemble approach ensures that 
these uncertainties are propagated into estimates of basin-mean and 
global-mean sea level. See ref. “ for further details about the GIA pre- 
dictions and the data used to weigh the GIA ensemble members. For 
each of our ensemble members, we randomly select one GIA prediction 
from the 128,000-member ensemble. Extended Data Fig. 4b shows 
the ensemble-mean RSL changes caused by GIA. Using the ICE6G D 
(VM5a) model* to account for GIA (Extended Data Fig. 5) does not cause 
noteworthy differences in global-mean and basin-mean observed sea 
level and the contributing processes. The differences in the subtropical 
North Atlantic basin are slightly larger (up to 0.3 mm yr”), but even here 
the GIA-related sea-level changes are within the confidence intervals 
of our GIA ensemble. 


Contemporary mass redistribution 

For the sea-level changes due to contemporary mass redistribution, 
we need to estimate the amount of water that is redistributed, and 
where on land the water is added or removed. During 2003-2018, we 
use GRACE and GRACE-FO observations, based on the JPL RLO6 mascon 
solution?’*“*, This solution provides monthly land-mass changes ona 
nominal 3-degree grid, from which we compute annual averages. Each 
grid cell has an associated measurement uncertainty, based on the for- 
mal error covariance matrix of the GRACE solution®. For each ensemble 
member, we randomly draw from these uncertainty estimates, perturb 
the mass estimates with this draw and correct for GIA. We then split the 
land-mass changes from GRACE into mass changes from glaciers, ice 
sheets and TWS using a previously described method*. 

Over 1900-2003, we use multiple estimates of each of the afore- 
mentioned processes. To combine these estimates with the GRACE 
observations, we average all observation-based mass-loss estimates 
over the same grid as the GRACE observations and remove the common 
mean in 2003 at every GRACE grid cell. Extended Data Fig. 6 shows all 
individual estimates and the resulting final composite estimate for 
each mass-redistribution process. 

For glaciers, we use two mass-change estimates. The first estimate, 
which covers the whole twentieth century, is based ona global glacier 
model thatis driven by observation-based surface forcing'®. This model 
produces estimates of the annual rate of glacier mass loss for each of 
the 19 glaciated regions defined in the Randolph Glacier Inventory 
(RGI)**. The second estimate”, which provides mass changes since 
1961, uses in situ glaciological and geodetic observations to derive 
total mass changes for each glaciated region. Both estimates provide 
uncertainties of the rate. For each ensemble member, we randomly 
choose between the two estimates. Before 1961, each member uses the 
estimates from the first estimate. Both estimates provide annual rate 


uncertainties. We draw random numbers using the rate uncertainties 
as the standard deviation, add them to the estimated rate and integrate 
this perturbed rate to obtain the total glacier mass changes. Because 
GRACE cannot distinguish the contributions from the Greenland and 
Antarctic peripheral glaciers from those from the ice sheets, we do not 
include these glaciers into the glacier mass balance. For Greenland, we 
add the peripheral glaciers to the ice-sheet contribution. For Antarctica, 
the mass balance ofits peripheral glaciers is very uncertain, owing to the 
lack of observations*’. However, since 2003, only a very small mass loss 
has been observed for these glaciers’, and observations since the 1950s 
donot suggest a large contribution”. Therefore, we assume no mass loss 
from the Antarctic peripheral glaciers. We account for missing (owing 
to their relatively small size) and disappeared glaciers using a previous 
estimate”. This study“ provides upper- and lower-bound estimates of 
the contribution of missing and disappeared glaciers. For each ensemble 
member, we uniformly sample between the upper- and lower-bound 
estimates. Since this estimate does not provide glacier mass changes per 
RGI region, we assume that the regional distribution of the contribution 
from missing and disappearing glaciers can be scaled by the regional 
relative contribution from the large glaciers as recognized by RGI. 

For the Greenland Ice Sheet, we use three estimates: a mass-balance 
reconstruction“ that covers 1900-2003, input-output estimates” 
that cover 1972-2003, and a multi-method assessment” that covers 
1993-2003. For each ensemble member, we randomly select one of 
these models. We use the first estimate for the contribution over the 
era for which the others do not provide an estimate. Each estimate 
provides a rate uncertainty, and we use these uncertainties to gener- 
ate a perturbed estimate for each ensemble member using the same 
procedureas for glaciers. These reconstructions (except the one from 
ref. *) do not include the contribution from peripheral glaciers. For 
these estimates, we add the estimated peripheral glacier contribution 
to the Greenland mass balance using the same approach as for other 
glaciated regions. 

For Antarctica, no mass-balance reconstruction exists before 
the satellite era, although observational evidence suggests 
twentieth-century mass loss, especially from West Antarctica. 
Therefore, we assume a small Antarctic Ice Sheet contribution before 
1993 of 0.05 + 0.04 mm yr”, based on an existing compilation”. For 
1993-2003, we use the multi-method assessments”*™ to derive the 
mass changes. To obtain an estimate of the spatial pattern of the mass 
changes from bothice sheets, we derive the spatial pattern of the mass 
loss from the perturbed GRACE observations. We assume this spatial 
pattern remains constant in time. 

The TWS component consists of natural and anthropogenic pro- 
cesses. For natural TWS, we use a twentieth-century reconstruction” 
that provides 100 ensemble members of natural TWS changes. We 
mask out all glacier and ice-sheet regions from these estimates, and 
randomly select one of the 100 TWS ensemble members. For anthropo- 
genic TWS changes, we consider artificial reservoir impoundment and 
groundwater depletion. For reservoir impoundment, we use an updated 
list of global artificial reservoirs”° and the ICOLDS dam database". We 
assume the filling and seepage rates of each reservoir follow previous 
estimates”. The ICOLDS dam database, which covers 93% of the total 
impounded volume, provides location coordinates of each reservoir; 
the database from ref.” does not. To approximate the regional distri- 
bution of this reservoir impoundment, we add the impounded water 
of the reservoirs with unknown location to the reservoirs with known 
location. We compute the fraction of the total impounded volume held 
by each known reservoir, and distribute the water from reservoirs with 
unknown location using this fraction. To our knowledge, for reservoir 
impoundment, no formal uncertainties have been quantified. Likely 
sources of the uncertainty in the reservoir impoundment stem from 
reservoir filling levels, storage-capacity loss due to sedimentation and 
seepage effects*”. Previous assessments assumed rates of uncertainties 
of 10%-30%*"8; we assume an uncertainty of 20% (10). 


For groundwater depletion, we use two gridded depletion estimates. 
Ref. °° provides depletion estimates over 1900-2003. However, a sub- 
stantial fraction of the depleted groundwater remains on land rather 
than ending up in the ocean”. To account for this effect, we assume 
that 40% of the depleted groundwater stays on land, and we scale the 
estimated depletion from this study by a factor of 0.6 (ref. **). We also 
use depletion estimates” over 1961-2003. Similarly to the glacier and 
ice-sheet case, we randomly select one of the estimates for each ensem- 
ble member. We assume an uncertainty of 20% (10) in groundwater 
depletion, which corresponds to previously estimated uncertainties”. 

These land-mass changes result in barystatic sea-level changes 
and, owing to GRD effects, in regionally varying sea-level change and 
solid-Earth deformation patterns. For each ensemble member, we 
solve the sea-level equation using a pseudo-spectral method>**. The 
spherical-harmonics transformations are computed using the SHTns 
library” up to degree and order 360. The resulting geoid changes 
and deformation are expressed relative to the centre-of-mass refer- 
ence frame, and include rotational feedback*®. We assume an elastic 
solid-Earth response to the land-mass changes, for which we use Love 
numbers based on® the Preliminary Referenced Earth Model®. With 
this procedure, we obtain 5,000 ensemble members, each consist- 
ing of annual time series of local sea-level changes and solid-Earth 
deformation due to contemporary mass redistribution. Extended Data 
Fig. 4a, c, e shows the ensemble-mean RSL trends due to contemporary 
mass redistribution. 


Steric changes 

We estimate global-mean and basin-mean steric changes for 1957-2018 
from gridded temperature and salinity reconstructions based onin situ 
observations of temperature and salinity. We use existing gridded esti- 
mates*’” for the upper 2,000 m. From these observations, we compute 
steric height anomalies using the TEOS-10 GSW software“. We also 
use gridded steric sea-level change estimates’. For each ensemble 
member, one of these estimates is selected randomly. Before the end 
of the 1950s, in situ observations are too sparse to derive unbiased 
steric changes”. For the upper-ocean (above 2,000 m) contribution 
before 1957, we use estimates computed from sea-surface temperature 
anomaly observations and estimates of ocean heat anomaly pathways 
from an ocean reanalysis. We also use the deep-ocean (below 2,000 m) 
steric expansion” for the full 1900-2018 period. These estimates come 
with an uncertainty, whichis used to perturb each ensemble member. 
In the Argo float data, a salinity drift has been detected since 20156, 
which causes an underestimation of global steric sea level. We correct 
for this drift by removing the estimated global-mean halosteric sea-level 
changes from each gridded estimate. Extended Data Fig. 7 depicts the 
time series of the individual steric products and the resulting estimates 
used in this paper. 


Sea-level observations 

We use annual-mean tide-gauge observations from the revised local 
reference (RLR) dataset from the Permanent Service for Mean Sea 
Level®**, as well as an extended tide-gauge dataset®, which has been 
updated until 2018. We remove observations that have been flagged 
for quality issues. Some stations show apparent data problems, such 
as spikes, jumps, drifts and large trends. These problems are typically 
caused by earthquakes, local subsidence, levelling issues and instru- 
ment problems. Owing to the multitude of the data problems, such 
stations cannot be automatically flagged and excluded, onthe basis of 
pre-set criteria, and we manually remove these regions from the analy- 
sis. We ultimately use 559 individual tide-gauge records in our recon- 
struction. From each sea-level record, we remove the self-consistent 
equilibrium nodal cycle” and the effects of local wind and sea-level 
pressure changes. To this end, we use wind and sea-level pressure fields 
from the ERA-20c reanalysis‘? over 1900-1979 and ERAS reanalysis® 
from 1979-2018, and use a simple linear regression model to remove 
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the wind and pressure effects”. Some locations, such as Aberdeen, 
Sydney and Singapore, have multiple tide-gauge records with dif- 
ferent observational periods. We merge stations that are within 
20 km of each other and have an overlap of at least 5 station years 
into regions. Henceforth, we refer to regions to denote any location 
that has a single or multiple merged tide-gauge observations. We 
only consider regions with at least 20 years of data. We link each 
region to a single ocean basin. All regions and the associated basins 
are shown in Extended Data Fig. 3. 


VLM 

Tide-gauge observations are affected by VLM”, and correcting these 
records for VLM has resulted in more coherent sea-level trends across 
different tide gauges”’’””. We use VLM observations from permanent 
GNSS stations and from the difference between satellite-altimetry 
and tide-gauge observations”. The RSL patterns associated with 
GIA and GRD are partially caused by solid-Earth deformation, which is 
observed as VLM. To avoid double-counting, we subtract the modelled 
solid-Earth deformation that results from GIA (Rca) and contemporary 
GRD (Rer) from the observed VLM time series (R,,,), to obtain a time 
series of residual VLM***: 


Reesidual(t) z Rops(t) g Royalt) = Rern(t) (1) 


We compute the linear trend in residual VLM, and we assume that 
the rate of residual VLM is representative for the full length of the 
tide-gauge record. 

We use the GNSS station database from the University of Nevada, 
Reno”. We select all GNSS stations that are within a 30-km radius of 
each region, have at least 4 years of daily observations, and for which 
the standard error of the residual VLM trend does not exceed 1mm yr”. 
We estimate the residual VLM trend using the MIDAS trend estimator”. 
We compute residual VLM for each ensemble member. The uncertainty 
in the derived trend is caused by the uncertainty in the corrections for 
GIA and contemporary GRD effects, and by the uncertainty that arises 
from estimating a linear trend from serially correlated data. The uncer- 
tainty due to GIA and contemporary GRD is estimated by computing the 
residual VLM trend for each individual ensemble member. To account 
for serial correlation, for each ensemble member we determine the 
trend uncertainty provided by the MIDAS trend estimator. We then 
draw a random number from a Gaussian distribution with this trend 
uncertainty as standard deviation, and perturb the estimated trend 
with this random number. 

To obtain residual VLM trends from the difference between 
satellite-altimetry and tide-gauge observations, we use the MEaSUREs 
gridded sea surface height anomalies version 1812 dataset”. This data- 
set has been corrected for calibration issues that caused a sea-level 
drift over the first years of the altimetry era**. The altimetry data covers 
the period 1993-2018. To obtain local residual VLM, we subtract GIA 
and contemporary GRD effects from altimetry. We require 15 years of 
overlap between altimetry and the tide gauge, and select all grid points 
within a300-km radius for which the correlation between annual-mean 
de-trended altimetry and tide-gauge sea level is above 0.5. This value, 
and the radius of 300 km, are chosen as a compromise between accu- 
racy and the number of locations for which VLM can be estimated”. 
We compute the residual VLM time series for each accepted altimetry 
grid point, and then compute the mean residual VLM time series by 
taking the mean of all individual time series, weighted by the correla- 
tion with the tide-gauge record. From this time series, we compute the 
linear trend and standard error by assuming that the serial correlation 
of the time series can be approximated by a first-order autoregressive 
process. This computation is performed using the Hector software”. 
For stations for which no single altimetry grid point has a correlation 
of 0.5 or higher, or for which the standard error is above 1 mm yr‘, no 
VLM estimate is generated. Similarly to the GNSS approach, we perturb 


each ensemble member with the trend uncertainty that arises from 
serial correlation in the time series. 

Some VLM observations appear as single outliers compared to nearby 
other observations, or result in unrealistically high or low sea-level 
trends. As for the tide-gauge selection procedure, owing to the mul- 
titude of possible problems in VLM estimates, no general criteria can 
be applied to catch these problems. Therefore, we manually remove 
VLM estimates that show such problems. For regions with multiple 
GNSS stations, or with both GNSS and altimetry VLM estimates avail- 
able, we use the average residual VLM trend, weighted by the inverse 
of the squared standard errors of the individual estimates. We are not 
able to estimate a VLM trend for all tide-gauge regions. For stations 
for which no VLM trend is available, we assume no residual VLM anda 
residual VLM standard error of 1 mm yr“. This standard error is based 
onthe maximum VLM uncertainty that we accept and on the stand- 
ard deviation among the residual VLM estimates, 1.5mm yr”. In some 
regions, large sea-level trends are compensated for by large residual 
VLMtrends. As a result, this standard deviation is probably biased high 
for regions without residual VLM estimates, because regions witha 
large sea-level trend and no residual VLM estimate are removed during 
the quality control phase. 


Global-mean and basin-mean sea-level reconstruction 

Following ref. °, before merging the individual region estimates into 
basin-mean curves, we estimate and remove the biases between local 
sea-level changes in each region and basin-mean sea-level changes that 
result from GIA, contemporary GRD effects and residual VLM. This 
correction results in an estimate of basin-mean sea level (Nasin), given 
observed regional sea level (7),egion), the difference between regional 
sea-level changes that result from GIA (fers region) ANd GRD (exp region), and 
the associated basin-mean sea-level changes, as well as residual VLM: 


Noasin)= N egion(É) T Isa, basin(É) z ncı, region(6)] 
t (erp, basin£) z Nero, region()] F R;esiduai(t) 


(2) 
Local sea-level variability may not be representative for the basin as a 
whole. To assess the uncertainty due to this non-representativeness, we 
perturb each ensemble member of the sea-level observations Nyegion(t) 
from each individual region with a realization of first-order autoregres- 
sive (ARI) noise. The ARI noise parameters are computed from the 
standard deviation and the first-order serial correlation of the regional 
sea-level observations. After computing all basin sea-level estimates 
from each individual region, we merge all the individual regions into 
asingle basin estimate using the virtual-station method®”°”’’, in which 
the two nearest regions are merged into a new virtual station halfway 
between the merged stations. Tide-gauge observations are not tied toa 
common vertical datum system. To account for different datum systems 
during the averaging process, we remove the common mean between 
two series estimated over their overlapping period. This procedure is 
repeated until one virtual station is left. The sea-level change estimate 
from the final virtual station is used as the basin-mean estimate. We 
obtain the final GMSL estimate by averaging the basin-mean estimates, 
weighted by the relative surface area of each basin. 

The resulting GMSL estimate shows a linear trend and multidecadal 
variability pattern that agree with other recent reconstructions*>””. 
These recent reconstructions all show lower twentieth-century rates 
than do earlier assessments’®*°, as shown in Extended Data Fig. 2. 

The global-mean and basin-mean altimetry curves are computed 
using the same gridded altimetry product as used for the VLM com- 
putations. To obtain basin-mean and global-mean RSL, we add the 
modelled deformation of the seafloor due to GIA and contemporary 
GRD effects to the altimetry curves®. 

The linear trends and accompanying uncertainty estimates in all 
basin-mean and global-mean quantities discussed here are computed 
from the linear trends in each ensemble member. Because the unique 


GIA model used in each ensemble member has an associated likelihood, 
we use the likelihood from the GIA model as the weight for the ensem- 
ble member when computing the mean and confidence intervals in all 
components. Because not all terms follow a Gaussian distribution, the 
confidence intervals are not assumed to be symmetric, and we directly 
compute the confidence intervals from the 5th and 95th percentile of 
the weighted ensemble. We account for the uncertainties due to serial 
correlation in the time series by adding the estimated trend uncertainty 
to the ensemble spread in quadrature. We assume that the spectrum 
of all time series can be approximated by a generalized Gauss-Markov 
spectrum™. We compute the noise parameters and the resulting trend 
uncertainty using the Hector software”. 


Data availability 


The resulting global and basin-scale reconstructions, the time series 
of global and basin sea-level changes and its contributors, grids with 
local sea-level and solid-Earth deformation due to contemporary GRD 
effects, and the individual ensemble members are available at https:// 
doi.org/10.5281/zenodo.3862995. 


Code availability 


The codes to compute the ensemble of observed sea-level changes and 
contributing processes, and the post-processing routines to compute 
statistics and to generate the figures are available at https://github. 
com/thomasfrederikse/sealevelbudget_20c. 
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Extended Data Table 1| Trends in observed basin-mean sea level and its contributors 


Subpolar North Atlantic 1900-2018 1957-2018 1993-2018 

Glaciers 0.42 0.32 0.52 0.31 0.21 0.41 0.36 0.28 0.47] 
Greenland Ice Sheet -0.16 -0.19 -0.13] -0.11 -0.14 -0.08] -0.25 -0.27 -0.22] 
Antarctic Ice Sheet 0.08 -0.01 0.17] 0.13 0.04 0.22 0.33 0.21 0.45] 
Terrestrial Water Storage -0.14 -0.23 -0.06] -0.08 -0.20 0.03] 0.16 0.04 0.28] 
Barystatic 0.19 0.05 0.34 0.25 0.09 0.41 0.61 0.40 0.82] 
Glacial Isostatic Adjustment 0.62 0.47 0.80 0.62 0.47 0.80 0.62 0.47 0.80] 
Steric - 0.62 0.40 0.86 1.18 1.00 1.37] 
Sum of Contributors - 1.50 1.10 1.92 2.42 2.00 2.82] 
Observed sea level 1.08 0.79 1.38 1.52 1.23 1.83 2.69 2.18 3.18] 
Altimetry - - 2.17 1.66 2.66] 
Indian Ocean-South Pacific 1900-2018 1957-2018 1993-2018 

Glaciers 0.73 0.55 0.90 0.56 0.35 0.72 0.73 0.56 0.89] 
Greenland Ice Sheet 0.48 0.39 0.58 0.33 0.24 0.42 0.72 0.63 0.80] 
Antarctic Ice Sheet 0.05 -0.02 0.12] 0.09 0.02 0.16 0.22 0.13 0.32] 
Terrestrial Water Storage -0.24 -0.37 -0.11] -0.17 -0.34 -0.01] 0.30 0.11 0.48] 
Barystatic 1.03 0.73 1.34 0.79 0.47 1.11 1.97 1.61 2.32] 
Glacial lsostatic Adjustment -0.15 -0.21 -0.08] -0.15 -0.21 -0.08] -0.15 -0.21 -0.08] 
Steric - 0.64 0.41 0.87 1.50 1.26 1.76] 
Sum of Contributors - 1.29 0.68 1.91 3.32 2.68 3.94] 
Observed sea level 1.33 0.80 1.86 1.51 1.03 2.00 3.93 3.32 4.55] 
Altimetry - - 3.65 3.23 4.08] 
Subtropical North Atlantic 1900-2018 1957-2018 1993-2018 

Glaciers 0.68 0.50 0.85 0.50 0.30 0.65 0.62 0.46 0.77] 
Greenland Ice Sheet 0.23 0.18 0.27 0.15 0.11 0.20 0.34 0.29 0.38] 
Antarctic Ice Sheet 0.10 -0.01 0.20] 0.16 0.06 0.26 0.40 0.26 0.52] 
Terrestrial Water Storage -0.13 -0.23 -0.03] -0.05 -0.18 0.09] 0.28 0.13 0.43] 
Barystatic 0.87 0.62 1.12 0.77 0.50 1.03 1.63 1.33 1.92] 
Glacial Isostatic Adjustment 0.76 0.40 1.04 0.76 0.40 1.04 0.76 0.40 1.04] 
Steric - 1.29 1.02 1.58 1.08 0.60 1.50] 
Sum of Contributors - 2.81 2.29 3.35 3.48 2.72 4.19] 
Observed sea level 2.49 1.89 3.06 2.76 2.05 3.42 3.98 2.75 5.20] 
Altimetry - = 4.04 2.77 5.24] 
East Pacific 1900-2018 1957-2018 1993-2018 

Glaciers 0.66 0.47 0.85 0.48 0.25 0.64 0.62 0.44 0.76] 
Greenland Ice Sheet 0.48 0.39 0.58 0.33 0.24 0.42 0.72 0.63 0.81] 
Antarctic Ice Sheet 0.09 -0.01 0.19] 0.15 0.05 0.24 0.37 0.24 0.49] 
Terrestrial Water Storage -0.22 -0.34 -0.09] -0.15 -0.32 0.02] 0.32 0.13 0.49] 
Barystatic 1.02 0.70 1.32 0.81 0.46 1.13 2.02 1.66 2.37] 
Glacial Isostatic Adjustment 0.03 -0.06 0.13] 0.03 -0.06 0.13] 0.03 -0.06 0.13] 
Steric = 0.47 0.21 0.74 0.37 -0.01 0.77] 
Sum of Contributors - 1.32 0.86 1.74 2.43 1.90 2.95] 
Observed sea level 1.20 0.76 1.62 1.64 1.26 2.03 1.82 1.10 2.56] 
Altimetry - - 2.35 0.70 4.06] 
South Atlantic 1900-2018 1957-2018 1993-2018 

Glaciers 0.76 0.56 0.95 0.56 0.32 0.73 0.72 0.54 0.88] 
Greenland Ice Sheet 0.50 0.41 0.60 0.34 0.24 0.43 0.74 0.65 0.83] 
Antarctic Ice Sheet 0.09 -0.01 0.19] 0.16 0.06 0.25 0.38 0.26 0.50] 
Terrestrial Water Storage -0.20 -0.35 -0.05] -0.12 -0.30 0.07] 0.38 0.18 0.58] 
Barystatic 1.15 0.82 1.48 0.93 0.58 1.26 2.23 1.87 2.61] 
Glacial Isostatic Adjustment -0.04 -0.08 -0.02] -0.04 -0.08 -0.02] -0.04 -0.08 -0.02] 
Steric - 0.88 0.73 1.03 1.29 0.98 1.66] 
Sum of Contributors - 1.78 1.23 2.33 3.48 2.79 4.15] 
Observed sea level 2.07 1.36 2.77 2.17 1.62 2.73 3.89 2.44 5.33] 
Altimetry - - 3.45 3.04 3.86] 
Northwest Pacific 1900-2018 1957-2018 1993-2018 

Glaciers 0.69 0.50 0.88 0.50 0.28 0.67. 0.65 0.47 0.80] 
Greenland Ice Sheet 0.53 0.42 0.64 0.35 0.26 0.45 0.78 0.68 0.88] 
Antarctic Ice Sheet 0.10 -0.01 0.20] 0.16 0.06 0.26 0.40 0.26 0.53] 
Terrestrial Water Storage -0.22 -0.36 -0.09] -0.16 -0.34 0.02] 0.34 0.15 0.52] 
Barystatic 1.09 0.76 1.41 0.86 0.51 1.21 2.17 1.78 2.53] 
Glacial Isostatic Adjustment -0.12 -0.22 -0.03] -0.12 -0.22 -0.03] -0.12 -0.22 -0.03] 
Steric - 0.71 0.45 0.98 1.20 0.79 1.58] 
Sum of Contributors - 1.44 0.80 2.07 3.23 2.52 3.93] 
Observed sea level 1.68 1.27 2.09 1.80 1.42 2.18 2.77 2.11 3.39] 
Altimetry - - 3.53 [2.64 4.45] 


The trends are given in millimetres per year, over 1900-2018, 1957-2018 and 1993-2018. The numbers in brackets indicate the 90% confidence intervals. 


